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EVAPORATION  INTO  COUETTE  FLOW 


1.  INTRODUCTION 

1.1  Objective 

The  objective  of  this  report  is  to  present  a  theory  for  the  evaporation  of  a 
small  drop  into  a  Couette  flow.  The  application  of  such  an  analysis  is  to  the 
experimental  data  obtained  from  the  Agent  Fate  5-cm  wind  tunnels  for  hazardous 
chemicals  and  comparison  with  some  of  that  data  is  provided  to  assess  the  usefulness 
of  the  analysis.  The  ultimate  purpose  being  to  develop  secondary  evaporation  models 
for  evaluating  the  effects  of  routine  industrial  emissions,  accidental  releases  of 
hazardous  materials,  and  dissemination  of  chemical  and  biological  warfare  agents. 

One  of  the  problems  in  modeling  of  the  evaporation  of  chemical  agents 
into  the  atmosphere  is  in  determining  the  way  evaporation  depends  on  the  wind  speed. 
Full-scale  atmospheric  tests  are  difficult  to  get  reliable  data  from  because  of  the  large 
number  of  variables  that  cannot  be  controlled.  Wind  tunnel  tests  can  provide  data 
under  controlled  conditions  and  theoretical  analysis  can  provide  the  relationships 
needed.  A  number  of  investigators  have  studied  these  problems  and  examining  some 
of  their  results  is  in  order. 

But  before  beginning  the  review  it  is  necessary  to  state  that  the  discussion 
will  be  limited  to  small  surface  drops  of  a  few  microliters.  Excluded  are  the  large-scale 
spills  where  normal  boundary  layer  techniques  are  more  applicable. 

Only  two-dimensional  flows  will  be  considered  here.  This  restriction 
means  that  the  evaporation,  subsequent  diffusion,  and  the  main  flow  velocity  vector 
remain  in  one  plane  normal  to  the  surface;  thus,  lateral  diffusion  and  cross  flows  are 
assumed  negligible. 

Here  the  primary  concerned  is  the  relationship  between  the  evaporation 
rate  and  the  variables  that  affect  it  such  as  drop  size  and  drop  properties.  Of  particular 
interest  is  the  effect  of  the  convective  external  flow  on  the  evaporation  rate.  Ultimately, 
there  is  also  concern  with  the  time-history  of  the  decreasing  mass  of  the  drop  and  its 
changing  geometry;  but,  this  is  not  addressed  here. 

1.2  Dimensional  Analysis 

Dimensional  analysis  can  be  used  to  present  the  evaporation  rate  in  terms 
of  three  dimensionless  parameters,  Sherwood  number  (Sh),  Reynolds  number  (Re), 
and  the  Schmidt  number,  (Sc).  The  set  of  basic  variables  are: 
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•  N  =  Evaporation  rate,  Kg/nrs 

•  N  =Average  drop  evaporation  rate,  Kg/nrs ,  =  M/A 

•  M  =  Total  drop  evaporating  mass  flux,  Kg/s 

•  A  =  Evaporating  surface  area,  m2 

•  D  =  Mass  diffusion  coefficient,  m2  /s 

•  v  =  Kinematic  viscosity,  m2  /s 

•  L  =  A  characteristic  length,  m 

•  u=  A  characteristic  velocity,  m/s 

•  cw  -cA  =  The  concentration  difference  in  Kg/m3  of  the  evaporating 

vapor  between  the  surface,  w,  and  the  edge  of  the  evaporating  vapor  dispersion  layer, 

A 

•  C  =  Proportionality  constant 

The  result  is: 

Sh  =  CRe"  Scm  (1) 


where 

N  •  L 

Sh  = - ,  Sherwood  number 

(cw  -ca)  D 


Re  =  — ,  Reynolds  number 

v 

Sc  =  — ,  Schmidt  number 
D 


1.3  Literature  Review 

Table  1  presents  the  results  of  a  brief  search  of  available  evaporation  rate 
estimates.  The  first  three  items  are  from  a  survey  by  Barry.1  It  is  interesting  that  none 
of  the  empirical  methods  are  concerned  with  the  diffusivity  of  the  evaporated  vapors. 
The  first  four  include  vapor  properties  through  the  surface  concentration  or  vapor 
pressure.  All  of  the  empirical  and  the  two  theoretical  boundary  layer  methods  apply  to 
large  spills.  A  problem  with  applying  the  large  spill  and  boundary  layer  approach  is  the 
fetch  or  boundary  layer  development  length  that  is  difficult  to  define  in  an  atmospheric 
environment.  Only  the  Baines  and  James2  and  the  present  prediction  specifically 
address  the  droplet  case. 
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The  Baines  and  James  paper  is  basically  identical  to  the  present  work, 
which  was  carried  out  independently.  The  only  differences  are  that  they  used  a 
similarity  solution  technique,  and  their  method  of  averaging  to  obtain  the  final  result  was 
to  match  the  circular  area  to  a  square  to  get  the  stream  wise  droplet  dimension  of 

L=J~z/4d . 


Table  1.  Evaporation  Rate  Formula  Cited  in  the  Literature 

Reference 

Evaporation  Rate 

n 

m 

Notes 

Empirical 

U  S  Air  Force*3 

Ncc  u0,5M.TF(Plp/P.p%) 

0.75 

Pvp  sat.  vapor 
pres. 

epa4& 

Nocu°'78M™67Pvp/RTA 

0.78 

Pvp  sat.  vapor 
pres. 

Stiver- 

MacKay67,8,9 

NccuP.(M„/RTa 

1.0 

P  sat.  vapor 

pres 

Sutton™ 

N  cc  u0'78L0'"(cw  -c„) 

0.78 

u=wind  vel. 
L=Free  liquid 
surface 

Coutant  and 
Penski11 

RcP=Rcp.0(l  +  cReHh/H)p) 

0.63 

Rcp  *  M 

u=Ave.  duct 
vel. 

H=duct  height 

Theoretical 

Laminar 

boundary  layer12 

ShL  =  0.66ReL^Sc^ 

1/2 

1/3 

u=free  stream 
L=b.l. 

development 

Turbulent 
boundary  layer12 

ShL  =0.036  ReL^  Sc^ 

4/5 

y2 

u=free  stream 
L=b.l. 

development 

Discontinuous 

Boundary 

conditions13 

Shd  =  0.334  Re^Sc^- 

"3f4Y,YdY/6" 

2 L 3 J  IlJ 

1/2 

1/3 

d=drop  dia. 
u=free  stream 

Vel. 

L=  b.l.  dev. 

Baines  &  James^ 

ShL  =0.840 

UxL 

V 

2 

3Sc1/3 

2/3 

1/3 

UT=Friction 

Vel. 

L=drop  length 

Present 

Shd  =  0.852 

utd 

V 

2 

3  Sc 1/3 

2/3 

1/3 

UT  =Friction 

Vel. 

d=drop  dia. 
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2. 


DROPLET  CONFIGURATION 


Figure  1  shows  measurements  of  three  droplets  on  glass.  The 
configuration  appears  to  closely  approximate  a  segment  of  a  sphere  depicted  as  the  arc 
of  a  circle.  In  free-fall  with  negligible  aerodynamic  effects,  the  drop  will  form  a  perfect 
sphere  as  its  minimum  energy  configuration.  In  fact  that  is  how  the  “initial”  diameters 
are  obtained. 


Typical  Drop  Conditions 
HD  on  Glass 


Initial  Drop  Diameter 


0.54  nun 


Final  Drop  Height 

Figure  1.  Configuration  of  Three  Sessile  Drops 

The  present  case  is  quite  different  in  that  the  drop  is  acted  upon  by  gravity 
on  a  solid  surface.  Surface  tension  must  still  be  a  part  of  picture.  A  new  effect  is  the 
action  of  surface  tension  in  the  vicinity  of  the  surface.  A  molecule  on  the  droplet  surface 
at  the  junction  of  the  liquid,  gas  and  solid  has  the  solid  molecule’s  attraction  to  contend 
with.  This  produces  a  contact-angle  that  may  be  considered  a  property  of  the  system. 


If  we  assume  the  shape  of  the  droplet  to  be  a  segment  of  a  sphere  then  the  radius  of 
curvature  is  dictated  by  the  contact  angle  and  the  volume  of  the  liquid.  Unfortunately, 
not  enough  is  known  about  the  properties  of  the  liquid  and  the  liquid  on  glass  to 
calculate  the  configuration  from  this  basic  approach  but  the  contact  angle  can  be 
calculated  from  the  given  measurements  of  surface  diameter  d  and  droplet  height  h.  All 
the  important  characteristics  of  the  droplet  can  be  defined  in  terms  of  d  and  the 

parameter  r|  =  — . 
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Figure  2.  Sketch  of  Droplet  Geometry 


First  the  radius  of  curvature  of  the  surface  droplet,  Rc ,  is 


(2) 


The  contact  angle,  p ,  is 


B  =  arctan  — 

Id 


-V)J 


(3) 


We  see  that  if  /?  is  a  constant  then  r|=  h Id  is  also  a  constant  independent 
of  the  volume.  The  average  of  the  experimental  data  is  r|  =  0.1046  ±  20%  and  the 

value  of  h/d  using  23.6°±4°  is  also  0.1046.  Using  this  value  for  h/d,  the  base  diameter 
(chord  in  cross  section)  of  the  spherical  segment  can  be  written  in  terms  of  h/d  and  the 
volume,  Q,  of  the  droplet. 


d  = 


Finally  the  surface  area  of  the  segment  can  be  calculated 


Aseg  =^-[4n2  +  l] 


(4) 


(5) 
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The  data  obtained  from  Figure  (1)  is  provided  in  Table  2. 


Table  2.  Measured  Values 

Q  cu  mm 

9 

6 

1 

Average 

d  mm 

6.29 

4.75 

2.74 

h  mm 

0.54 

0.58 

0.29 

h/d 

0.0858 

0.1221 

0.1058 

0.1046 

/?  Deg.  Equation 
(2) 

19.5 

27.4 

23.0 

23.6 

R  mm  Equation  (1) 

9.428 

5.153 

3.381 

The  calculations,  based  on/?,  for  the  three  droplets  are  given  in  Table  3. 

Included  are  calculations  of  the  wetted  area,  7id2/4,  which  shows  the  relatively  small 
error  that  would  be  made  in  neglecting  the  curvature  of  these  small  droplets.  This  also 
provides  a  check  that  the  curved  surface  area  is  correctly  calculated. 


Table  3.  Calculated  Values 

Q  cu  mm 

9 

6 

1 

Average 

d  mm  Equation  (4) 

5.999 

5.2411 

2.884 

d  (meas.)/d  (cal) 

1 .0484 

0.9063 

0.8500 

0.9682 

Rc  mm,  Equation  (1) 

9.428 

5.153 

3.381 

A^sq  mm, 

Equation  (5) 

30.32 

23.14 

7.01 

A  Droplet  Wetted 
Area,  sq  mm 

28.27 

21.57 

6.53 

Since  the  diameter  in  the  segment  area  formula  is  proportional  to  a 
segment  volume  raised  to  the  1/3  power  then  the  area  of  the  segment  is  proportional  to 
the  2/3  power  of  the  volume.  The  proportionality  factor  is  just  a  function  of  h/d.  We  can 
write  the  area  as 


A  =f  — 


seg 


(6) 


The  result  is  plotted  in  the  Figure  3.  If  the  contact  angle  varied  by  5 
degrees,  the  f(r|)  could  vary  by  20%  and  the  evaporation  rate  changes  proportionatly. 
Data  on  the  effect  of  temperature  on  contact  angle  of  HD  are  not  currently  available. 
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0  10  20  30  40  50  60  70  80  90 

Contact  Angle,  Deg, 

Figure  3.  Shape  Factor  as  a  Function  of  Contact  Angle 


3.  DERIVATION  OF  INTEGRAL  THEORY 

In  this  derivation  the  protuberance  of  the  droplet  into  the  boundary  layer 
will  be  neglected.  The  largest  drop  in  the  thinnest  boundary  layer  considered  for  the 
5-cm  tunnel  tests  extends  as  far  into  the  flow  as  the  thickness  of  the  linear  region  of  the 
sublayer.  Because  the  droplet  can  be  considered  slender  (d/h  =  10)  the  boundary  layer 
can  ride  smoothly  over  the  drop.  However,  the  magnitude  of  the  approximation  made 
by  neglecting  the  disturbance  has  not  been  determined. 

It  is  also  assumed  that  the  evaporation  rate  is  determined  by  the 
concentration  distribution  in  the  vicinity  of  the  droplet  and  that  the  vapors  from  the 
droplet  are  confined  to  the  linear  sublayer  at  least  as  an  approximation.  It  is  further 
assumed  that  the  vapors  added  are  so  dilute  that  they  do  not  change  the  composition  of 
the  main  gas  flow  or  produce  a  significant  change  in  the  sublayer  momentum  flow  field. 

For  example  the  displacement  effect  of  the  added  gas  does  not  change 
the  velocity  field.  The  sublayer  is  modeled  as  Couette  flow  with  the  same  normal 
velocity  gradient  as  the  turbulent  sublayer. 


( du  ' 

2 

(u  } 

UH 

ldy.; 

V 

Turb,sl 

l  H  J 

Couette 


(7) 
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where  uH  is  the  moving  upper  wall  velocity  at  a  distance  H  from  the  fixed  wall.  The 
same  integral  procedure  as  that  used  for  the  temperature  step  boundary  condition 
problem  in  heat  transfer12  is  used  here. 


Uh 


*■ 


Figure  4.  Sketch  of  Droplet  in  Couette  Flow 


The  first  element  is  to  define  the  concentration  thickness,  3. 


A 

(cw  -cH)uHd=  Jucdy 
o 


(8) 


It  should  be  noted  that  cH  just  denotes  the  concentration  of  the 
evaporated  vapor  present  in  the  ambient  air  but  not  coming  from  the  droplet  (thus,  in 
this  case,  it  is  zero).  It  is  assumed  uniform  in  all  the  oncoming  flow. 

The  velocity  distribution  is 


u  =  uH 


_y 

H 


(9) 


In  this  boundary  analysis  the  droplet  and  the  concentration  distribution  are 
taken  as  two-dimensional  to  simplify  the  development.  The  concentration  distribution  is 
assumed  as  in  the  simple  laminar  boundary  layer: 


c  =  (cw  -cH)  1 


Ay+ilY| 

2  A  2{a) 


(10) 


14 


Let  <;  =  —  where  A  is  the  penetration  depth  of  the  vapor  concentration.  Combine  these 

A 

relations  so  that  9  becomes 


The  next  step  is  to  relate  the  rate  of  change  of  9  with  respect  to  the 
longitudinal  coordinate  x  with  the  wall  boundary  condition.  Note  it  is  assumed  that  x  =  0 
defines  the  start  of  the  two-dimensional  drop. 


(cw  -CH )uh  T“  = 
dx 


dc  ^ 


dy>U 


-  N 


(12) 


where  D  is  the  diffusion  coefficient  of  the  vapor  in  air  and  N  is  the  mass  flux  per  unit 
area  of  vapor  from  the  fixed  wall.  Using  the  definitions  of  9  and  c  and  recognizing  that 
A(x)  is  the  only  function  of  x. 


I_U]La— 

5  H  dx  “  2  A 


(13) 


This  differential  equation  can  be  solved  using  the  initial  condition  that  A  =0  at  x=0  to 
give 


Figure  5.  Sketch  Depicting  the  Flow  of  Evaporated  Vapor 

The  next  step  is  to  use  this  result  (Equation  (14))  in  the  definition  of  9 
(Equation  (11))  and  put  these  into  the  left-hand-side  of  Equation  (12)  to  define  N  . 
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N  =  C(cw  -ch)D 


=  0.531 


(15) 


Dx 

-X  2 

•  c  -  — 

f45l 

1 

c 

33 

i _ 

30 

l  2  J 

Equation  15  can  be  put  into  nondimensional  form  of  an  Sherwood  number, 
Sh  d ,  where  the  length  scale  is  the  diameter  of  the  drop,  d .  The  right-hand-side 

contains  a  sublayer  Re  number  where  we  will  use  Equation  (7)  to  replace  uh/H  by 
uj/v .  The  diffusivity,  D,  divided  into  the  air  kinematic  viscosity  forms  the  Sc  number, 

Sc  =  v/D . 


Shd  = 


Nd 


=  C 


u,d 

V 

Id  J 

-1/3 


(16) 


(Cw  "Ch)D 

Evaluation  of  an  average  Sh  may  be  obtained  by  integrating  over  x/d  from 


Oto  1. 


2 

3  Sc1 1/3 ,  |c  =  0.796  (17) 

This  is  the  simple  averaging  over  the  droplet  diameter  or  taking  the 
averaging  on  the  centerline  of  the  actual  drop.  The  actual  drop  is  circular  in  wetted  area 
on  the  surface. 


Shd  = 


Nd  _  3  c 
(cw-ch)D  2 


M 
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A  better  approximation  to  the  effect  of  the  circular  geometry  is  to  consider 
integrating  over  the  actual  area  of  the  drop  as  illustrated  in  Figure  6. 


7td2N 


/ 1 

=  M  oc  2  JjNdxdz 


(18) 


0  0 


the  z  coordinate  is  in  the  surface  perpendicular  to  the  normal  x-y  plane. 


where 

And  if  N  oc 


2 

l 

fdl 

,2, 

+  z“  = 

v2> 

(19) 


\  /3 

X  ' 


U  J 


then  N  =  1 .07  N  or  in  nondimensional  terms. 


Shd  =0.852 


Sc 


1/3 


(20) 


This  is  the  final  result,  which  shows  the  appropriate  Re  number  can  be 
based  on  the  shear  velocity  and  the  droplet  diameter. 

The  ratio  of  A  at  the  trailing  edge  of  the  drop  to  the  actual  thickness  of  the 
laminar  sublayer  can  now  be  computed.  The  linear  sublayer  thickness  is  taken 
asy+  =  4 ,  corresponding  to  a  1  %  deviation  from  the  linear  velocity  distribution.  If  we 
call  this  point  5S|  then 


f45l 

A 

Sc~1/3 

1 

O 

\  ir«i 

l  2  J 

V 

4  V  2  J 

0.706 


(21) 


4.  CFD  COMPUTATION  OF  COUETTE  FLOW  EVAPORATION 

The  numerical  computation  of  the  evaporation  from  a  sessile  drop  into  a 
Couette  flow  has  been  undertaken  as  a  first  step  in  applying  numerical  techniques  to 
the  case  of  turbulent  boundary  layers  including  three-dimensional  diffusion  effects.  This 
problem  is  basically  the  same  as  described  for  the  integral  method. 

The  numerical  code  used  is  a  modification  of  a  program  published  by 
Morgan14  for  computation  of  unsteady  one-dimensional  heat  transfer.  The  solution 
technique  is  an  implicit  Crank-Nicolson  method  with  two-point  boundary  condition  and 
marches  in  x  from  an  initial  condition. 
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4.1 


Diffusion  Equation 


The  equation  to  be  studied  here  is  the  diffusion  equation,  which  for  the 
two-dimensional  case  in  a  Couette  flow  is: 

3c  3~c 

u —  =  D— ,  on  the  drop  c(x,y=0)=cw  ,  c(x,  y  go  )=0  (22) 

dx  dy 

where  y  is  the  distance  normal  to  the  surface  and  x  is  the  distance  down  stream  on  the 
drop.  D  is  the  diffusivity  and  u  is  assumed  linear  in  y  for  Couette  flow. 

u  - 

In  such  a  flow  u  =  —  y ;  thus,  the  governing  equation  differs  from  the  heat 

v 

transfer  equation  because  of  a  variable  coefficient  involving  y.  Since  the  penetration 
depth  of  the  vapor  is  in  practice  finite,  the  second  boundary  condition  is  simplified  to  a 
sufficiently  large  constant  y. 

Equation  (22)  can  be  put  into  the  form: 


dC  d2C 
Y  dX  "  dY2  ’ 


(23) 


where  Y^,  X=^^,andC=^- 
10v  lOOOvSc  cw 

The  boundary  conditions  are  Y=0,  C=1  at  the  surface  of  the  drop  and  Y  = 

1 ,  C=0  at  y  =10  assuming  this  value  of  y  is  sufficient  away  from  the  vapor  plume. 
The  initial  condition  at  X  =  0,  C=0  for  all  Y  except  at  Y=0,  C=1  corresponding  to  a  step 
in  the  initial  boundary  condition  at  the  drop. 

4.2  Finite-Difference 

Equation  (23)  can  be  written  in  finite-difference  form  using  central 
differences  in  both  directions  but  with  the  node  taken  at  the  mid  point  of  the  X  step. 
Taking  i  to  be  the  ith  step  in  X  and  j  in  Y. 


fc  ,  -C 

1+1,  j  1,J 

_  i 

( r  -2  C  +  C  ^ 

^i+l,j+l  ^^i+lj  '  ^i+lj-l 

1 

J  

fCi,J.,-2C„J+C,H'| 

> 

X 

1 

1 

(AY)1  J 

n 

2 

l  (AY)2  J 

The  unknowns  are  in  the  i+1  X  step  where  there  are  three  unknown 
Cj+1  j_, ,  Ci+1J ,  and  Cj+I  j+l .  The  square  bracketed  term  is  the  only  addition  to  the  code 

developed  by  Morgan.  It  is  convenient  to  define  a  combined  step  variable  r. 
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AX 


(25) 


(AY)2 

The  three  unknowns  have  coefficients  a,  b,  c  and  the  remainder  d. 


a  =  -— ,  b  =  Yj  +r ,  and  c  =  (26) 

2  2 

d  =  (Y1-r)C11+i(CiJ.l+Ci.H)  (27) 


A  tridiagonal  matrix  is  constructed  from  the  a,  b,  c  and  d’s  for  j=2  to  N-1 
and  the  two  boundary  conditions  at  j=1  and  N.  The  Thomas  algorithm  can  be  used  to 
determine  all  of  the  N-2  unknowns  in  C  at  each  X  step  and  as  the  solution  marches 
downstream  the  entire  concentration  distribution  is  defined. 


4.3  Stability 

The  Crank-Nicolson  implicit  scheme  is  stable  for  any  step  size  in  a  wide 
range  of  problems  that  vary  smoothly  and  not  too  rapidly.  It  may  be  viewed  as  a  Fourier 
series  representation  of  the  solution  where  if  the  function  varies  rapidly  or  is 
discontinuous.  The  series  involves  many  high  frequencies.  It  is  the  high  frequencies 
that  are  unstable  in  the  Crank-Nicolson  method.  A  way  of  making  it  more  stable  is  to 

reduce  the  step  to  meet  the  stability  criterion  of  the  explicit  techniques  that  is  r  -  0  5 .  As 
the  solution  moves  away  from  a  discontinuity  the  step  size  can  be  increased  to  a  more 
efficient  value. 


4.4  CFD  Results 

The  application  of  the  Morgan  code  is  to  that  of  a  two-dimensional  drop  of 
HD,  which  extends  0.006  m  stream-wise  in  a  Couette  flow  of  friction  velocity, 
ut=0.2  m/s  and  kinematic  viscosity  of  0.0000159  m2  /s .  The  surface  concentration  of 

HD  is  0.002  Kg/ m3 and  Sc  number  is  2.53  from  which  the  diffusivity  of  HD  in  air  can  be 
calculated. 


The  drop  stream-wise  dimension  (designated  d  in  analogy  with  the  drop 
surface  diameter)  is  considered  the  maximum  distance  in  x  to  be  computed.  When  it  is 
put  into  the  nondimensional  form  suggested  by  Equation  (23),  X  is  0.02983  at  d.  Initially 
40  steps  were  selected,  which  results  in  an  AX=0.00075.  A  y+  of  10  was  thought  to  be 
sufficient  to  define  the  outer  boundary  condition  of  c=0.  Two  hundred  steps  in  Y  were 
assumed  or  AY  =  0.005. 
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The  evaporation  rate,  N  (in  Kg/m2s),  can  be  determined  by  Fick’s  Law. 


N  =  -D 

V 


dc 

dy 


\ 

/ 


(28) 


In  terms  of  the  computational  variables,  Equation  (28)  can  be  most  easily 
rewritten  in  nondimensional  Sherwood  number  terms.  (Note  the  length  d  makes 
Shd  nondimensional  but  N  and  Shd  are  still  functions  of  x.) 

Shd  =  =  Red[C(l)  -  C(NYOUT)]/[10(NYOUT  -  1)DY]  (29) 

c„,D 


where  Re  d  =  u  td  /  v  and  NYOUT  is  the  j  step  used  to  calculate  the  concentration 
derivative  (>2  and  is  not  sensitive  to  the  choice  because  of  the  linearity  near  the 
surface.)  Figure  7  illustrates  the  distribution  of  Shd  versus  the  computational  X.  The 

line  drawn  through  the  data  was  calculated  by  taking  a  point  near  X=0.03  and  running  a 

_  1/ 

curve  through  that  point  proportional  to  X  /3  as  predicted  by  the  integral  technique  and 
the  similarity  theory  of  Baines  and  James.13 


0.00  0.01  0.01  0.02  0.02  0.03  0.03  0.04 

x 

Figure  7.  Local  Sherwood  Number  versus  Computational  X 


Figure  8  more  definitively  confirms  that  the  CFD  results  duplicate  the 
results  of  the  integral  analysis  and  that  of  reference  2.  The  line  designated  as  “Integral” 
is  Equation  (16).  The  CFD  data  differ  from  Equation  (16)  only  at  the  smallest  values 


20 


of  x/d,  but  this  is  the  region  of  the  steepest  gradient  in  the  distribution  and  the  region  of 
the  greatest  uncertainty  in  the  numerical  computation. 


Concentration  profiles  can  be  obtained  from  the  numerical  computation 
and  examples  of  a  few  profiles  are  given  in  Figure  9.  The  four  profiles  are  specified  in 
terms  of  the  numerical  x  variable  X  but  the  corresponding  x/d’s,  i.e.,  the  distance  from 
the  leading  edge  in  percent  of  the  assumed  drop  diameter  (6  mm)  are  given  in 
Table  4. 


Table  4.  Comparison  of  CFD  and  Integral  Theories 

X 

x/d  (%) 

7  (%) 

a 

Shd  (CFD) 

Shd 

(Integral) 

0.000427 

1.43 

2.8 

53.2 

53.2 

0.00157 

5.25 

4.3 

34.5 

34.5 

0.00817 

27.40 

7.5 

19.9 

19.9 

0.0307 

100.00 

11.6 

12.8 

12.9 

Also  shown  on  Figure  9  are  profiles  computed  from  Equation  10  of  the 
integral  theory.  And  the  penetration  thickness  A  is  given  by: 
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A  =  d 


45  1  x 


(30) 


2  Sc  Re;;  d 


which  is  a  rewritten  form  of  Equation  (14).  Values  of  the  penetration  thickness  as  a 
percent  of  the  diameter  are  tabulated  in  Table  4. 


Figure  9.  Concentration  Profiles 


It  is  apparent  from  Figure  (9)  that  the  integral  assumption  closely 
approximates  the  correct  CFD  profile.  This  is  confirmed  by  comparison  of  the 
calculated  Sh  numbers  that  may  be  viewed  as  nondimensional  normal  concentration 
gradients  at  the  surface  (see  Equations  (28)  and  (29).)  The  calculated  Sh  numbers 
from  the  CFD  and  Integral  theories  are  given  in  Table  4  and  are  in  almost  perfect 
agreement. 


5.  COMPARISON  WITH  EXPERIMENTAL  DATA 

Information  is  available  which  provides  data  on  HD  diffusivity  and  its 
partial  pressure  as  a  function  of  temperature.15 

5.1  Properties  of  HD 

The  values  determined  are  given  in  Table  5.  To  convert  partial  pressure 
into  concentration,  the  molecular  weight  is  required  and  a  value  based  on  the  chemical 


22 


formula  C4H8C12S,  which  gives  a  Molecular  weight  =  159.07  Kg/Kmole.  Thus,  the  gas 

Nm 

constant  becomes  11®=  52.27  - .  The  perfect  gas  formula,  Equation  (31),  was 

KgK 

used  to  calculate  the  concentration  listed  in  Table  5. 


Phd 

R-hdT 


(31) 


where  the  pressure,  pHD ,  is  the  saturation  vapor  pressure  of  HD.  The  temperature  is 
the  absolute  ambient  temperature.  The  concentration  of  HD  in  the  free-stream  is  taken 
as  zero.  Also  given  in  Table  5  are  the  values  of  Sc  number. 


Table  5.  Properties  of  HD 

Temperature 

Vapor 

Pressure 

Concentration 

Diffusivity 

Sc  No. 

°C 

N/m2 

Chd,  Kg/m2 

m2/s 

Sc 

0 

1.32 

15 

5.01 

0.000333 

5.50E-06 

2.65 

20 

9.20 

25 

14.66 

35 

32.4 

0.002013 

6.50E-06 

2.53 

40 

46.66 

50 

105.4 

0.00643 

7.00E-06 

2.55 

5.2  Comparison  with  Data 

Measurements16  are  available  of  the  total  evaporation,  M ,  in  micro  g/min 
in  the  5x5  cm  tunnels  versus  the  nominal  free-stream  velocity  for  three  droplet  volumes 
(9,  6,  1  cu  mm)  and  for  three  tunnel  system  temperatures  (15,  35,  50  °C.)  The  drops 
considered  here  are  all  for  HD  on  a  glass  surface.  Tabulated  data  are  provided  in  the 
attached  Table  6.  The  areas  and  diameters  of  the  droplet  are  taken  from  previous 
discussion  of  the  droplet  configuration.  Sutherland’s  dynamic  viscosity  formula  was 
used.  The  density  was  calculated  from  the  specified  temperature  and  standard 
atmospheric  pressure  of  101 .325  KPa.  In  Table  5  the  average  friction  velocity  for  the 
three  speed  ranges  were  taken  from  the  recent  evaluation  of  the  velocity  profiles17  in  the 
5-cm  tunnels. 
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The  Sh  number,  Sh ,  data  are  graphed  versus  the  combined  parameter, 

Re^Sc^,  in  Figure  10,  with  temperature  as  a  parameter,  including  the  over-all 
correlation  line  (linear  statistical  best  fit.) 

The  slope  b  of  Sh  (=yj)  versus  Re 2/3  Sc1/3(=xi)  has  been  computed  using 
the  conventional  method  for  a  correlation  assumed  to  go  through  the  origin:  that  is 
summing  over  the  N=19  data  pairs: 

lyi 


i 


and  the  standard  deviation  in  the  y  data  relative  to  the  best  fit  line  is: 
<T,=^Z(  yj-bx;)2  (33) 


The  standard  deviation  in  the  slope  is  obtained  by  the  propagation  of  error 

method. 


=  G, 


(34) 
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The  experimental  result  is  a  slope  of  0.98  ±17%;  whereas,  Equation  (20) 
predicts  that  the  slope  of  the  line  should  be  0.852  or  13%  smaller.  However,  it  should 
be  kept  in  mind  that  the  analysis  is  for  the  simpler  two-dimensional  case  whereas  the 
drop  concentration  distribution  is  actually  three-dimensional.  That  is  the  concentration, 
in  the  three  dimensional  case,  is  free  to  defuse  latterly  from  the  drop  as  well  as  normal 
to  the  wall. 


Also  on  the  figure  is  the  uncertainty  band  based  on  plus  or  minus  1 .96 
sigma  or  a  confidence  level  of  95%.  The  upper  and  lower  slope  of  this  band  are  ±  17% 
of  the  best  fit  value  of  0.98  and  is  affected  by  some  of  15  °C  data,  which  fall 
considerably  outside  of  the  band.  If  these  points  were  omitted  from  the  statistical 
analysis  the  resulting  best-fit  slope  would  be  0.83  or  about  2.5%  lower  than  the 
predicted  value. 

Some  of  the  inaccuracies  may  be  due  to  the  vapor  pressure  and  diffusivity 
values  that  were  obtained  from  graphical  data.  Or  the  uncertainty  in  the  actual  drops 
size  and  area  that  may  be  affected  by  temperature  that  has  not  been  accounted  for  in 
these  results.  There  is  scatter  in  all  the  measurements  but  somewhat  more  in  the  15  °C 
data  as  might  be  expected. 


6.  CONCLUSIONS 

By  developing  a  theory  for  the  evaporation  from  a  droplet  under  conditions 
appropriate  to  the  Fate  wind  tunnel  an  attempt  has  been  made  to  show  that  a  relatively 
simple  Couette  diffusion-convection  model  could  provide  a  correlation  formula  that 
would  generalize  the  specific  test  measurements  (at  different  droplet  volume, 
temperature  and  different  sublayer  shear  velocity.) 

The  CFD  code,  modified  to  solve  the  2-D  Couette  flow  evaporation 
problem,  provides  identical  results  as  the  integral  theory  and  the  similarity  theory  of 
Baines  and  James. 

The  CFD  calculation  provides  reliable  details  to  the  solution  unavailable  to 
the  approximate  theories 

The  CFD  approach  has  the  prospect  of  extension  to  three  dimensions  and 
into  a  turbulent  boundary  layer. 

The  resulting  correlation  as  exhibited  in  Table  6  and  Figure  10  bring  the 
Couette  theory  into  a  certain  level  of  agreement  with  the  experimental  data.  The  slope 
in  the  statistical  best-fit  line  is  0.98,  which  is  13%  higher  than  the  Couette  theory’s  value 
of  0.852.  The  theory’s  slope  falls  slightly  above  the  lower  2-sigma  uncertainty  band  and 
if  the  three  high  15  °C  data  were  omitted  as  outliers  the  best-fit  slope  would  be  0.83  or 
2.5%  below  the  prediction. 
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There  seems  to  be  a  slight  trend  in  the  data  at  the  higher  velocities:  that 
is,  higher  Re  numbers  tend  to  increase  faster  than  the  linear  correlation. 

This  model  indicates  that  the  shear  velocity  and  droplet  diameter  are  the 
appropriate  length  and  velocity  scales  for  the  evaporation  in  this  case.  There  is  also 
reason  to  think  that  the  temperature  effects  can  be  understood  through  the 
incorporation  of  the  temperature  dependent  thermodynamic  properties  of  the  droplet 
material. 


Table  6a.  Input  Data 

T 

Q 

M  dot 

u  delt 

u  tau 

A 

d 

D 

nu 

c  surface 

C 

milli  1 

micro 

g/min 

m/s 

m/s 

sq  mm 

m 

sq  m/s 

sq  m/s 

Kg/cu  m 

15 

9 

4.1 

0.22 

0.0393 

30.3 

0.006 

5.50E-06 

1 .46036E-05 

0.000333 

15 

9 

19 

3.61 

0.181 

30.3 

0.006 

5.50E-06 

1.46036E-05 

0.000333 

15 

6 

3.1 

0.2 

0.0393 

23.1 

0.00524 

5.50E-06 

1.46036E-05 

0.000333 

15 

6 

11.7 

1.61 

0.096 

23.1 

0.00524 

5.50E-06 

1 .46036E-05 

0.000333 

15 

1 

1.5 

0.26 

0.0393 

7.01 

0.00288 

5.50E-06 

1 .46036E-05 

0.000333 

15 

1 

5.6 

3.61 

0.181 

7.01 

0.00288 

5.50E-06 

1 .46036E-05 

0.000333 

35 

9 

53.7 

1.6 

0.096 

30.3 

0.006 

6.50E-06 

1 .64466E-05 

0.002013 

35 

6 

20 

0.2 

0.0393 

23.1 

0.00524 

6.50E-06 

1 .64466E-05 

0.002013 

35 

6 

37.6 

1.61 

0.096 

23.1 

0.00524 

6.50E-06 

1 .64466E-05 

0.002013 

35 

6 

63.2 

3.58 

0.181 

23.1 

0.00524 

6.50E-06 

1 .64466E-05 

0.002013 

35 

1 

14.5 

1.61 

0.096 

7.01 

0.00288 

6.50E-06 

1 .64466E-05 

0.002013 

50 

9 

51.2 

0.22 

0.0393 

30.3 

0.006 

7.00E-06 

1.78814E-05 

0.00643 

50 

9 

196 

1.6 

0.096 

30.3 

0.006 

7.00E-06 

1.78814E-05 

0.00643 

50 

9 

272 

3.61 

0.181 

30.3 

0.006 

7.00E-06 

1.78814E-05 

0.00643 

50 

6 

42.8 

0.2 

0.0393 

23.1 

0.00524 

7.00E-06 

1.78814E-05 

0.00643 

50 

6 

100.2 

1.61 

0.096 

23.1 

0.00524 

7.00E-06 

1.78814E-05 

0.00643 

50 

1 

21.6 

0.26 

0.0393 

7.01 

0.00288 

7.00E-06 

1.78814E-05 

0.00643 

50 

1 

45.6 

1.61 

0.096 

7.01 

0.00288 

7.00E-06 

1.78814E-05 

0.00643 

50 

1 

66.9 

3.61 

0.181 

7.01 

0.00288 

7.00E-06 

1.78814E-05 

0.00643 
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Table  6b.  Nondimensional  Parameters 

T 

Q 

Sh= 

Re 

ReA2/3 

Sc 

ScA1/3 

R2/3S1/3 

(yi-b*xi)A2 

C 

milli  1 

yi 

xi 

15 

9 

7.39E+00 

16.14666 

6.388347 

2.66E+00 

1 .384733 

8.846156 

0.391307 

15 

9 

3.42E+01 

74.36503 

17.68395 

2.66E+00 

1.384733 

24.48755 

145.3124 

15 

6 

6.40E+00 

14.10142 

5.836805 

2.66E+00 

1.384733 

8.082418 

0.85128 

15 

6 

2.42E+01 

34.44621 

10.58671 

2.66E+00 

1.384733 

14.65977 

118.189 

15 

1 

5.61  E+00 

7.750397 

3.91636 

2.66E+00 

1 .384733 

5.423114 

0.483348 

15 

1 

2.09E+01 

35.69521 

10.8411 

2.66E+00 

1.384733 

15.01203 

53.83475 

35 

9 

1.35E+01 

35.02249 

10.70446 

2.53E+00 

1 .36266 

14.58653 

0.109575 

35 

6 

5.78E+00 

12.52127 

5.392195 

2.53E+00 

1 .36266 

7.347726 

0.76985 

35 

6 

1.09E+01 

30.58631 

9.780282 

2.53E+00 

1 .36266 

13.32719 

1.461168 

35 

6 

1.83E+01 

57.66794 

14.92644 

2.53E+00 

1  36266 

20.33966 

0.027031 

35 

1 

7.59E+00 

16.8108 

6.562342 

2.53E+00 

1 .36266 

8.942238 

0.262752 

50 

9 

3.75E+00 

13.18687 

5.581631 

2.55E+00 

1.366998 

7.630081 

9.971887 

50 

9 

1.44E+01 

32.2122 

10.12388 

2.55E+00 

1.366998 

13.83933 

3.365859 

50 

9 

1.99E+01 

60.73341 

15.45083 

2.55E+00 

1.366998 

21.12126 

0.657115 

50 

6 

3.60E+00 

11.51653 

5.099737 

2.55E+00 

1.366998 

6.971333 

7.399759 

50 

6 

8.42E+00 

28.13198 

9.249826 

2.55E+00 

1.366998 

12.6445 

9.230393 

50 

1 

3.29E+00 

6.329697 

3.421805 

2.55E+00 

1.366998 

4.677602 

0.905171 

50 

1 

6.94E+00 

15.46185 

6.206418 

2.55E+00 

1.366998 

8.484165 

0.560445 

50 

1 

1.02E+01 

29.15204 

9.472093 

2.55E+00 

1.366998 

12.94834 

2.409719 

sum  yi= 

2.25E+02 

sigmaY= 

4.448426 

sum  xi= 

229.371 

b(slope)= 

9.82E-01 

sigma  b= 

0.082282 
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NOMENCLATURE 


A  Surface  Area  of  evaporating  drop  (also  used  to  denote  an  matrix) 

a,  b,  c,  d  Coefficients  in  the  tridiagonal  matrix 
b  Best-fit  slope 

C  Proportionality  constant  or  Computer  variable=c/cw 

cH  Ambient  concentration  not  from  the  drop  (assumed  zero  in  this  study) 

cw  Concentration  of  vapor  at  the  surface  of  drop 

d  Diameter  of  sessile  drop 

D  Diffusivity  of  evaporated  vapor  in  air 

H  Channel  height 

h  Height  of  sessile  drop 

i  x-direction  chord  of  circular  surface  droplet  at  z  from  centerline 

L  Characteristic  length 

M  Total  drop  evaporation  rate 

Mw  Molecular  weight 

n  Number  of  data  points 

N  Local  evaporation  rate  per  unit  area 

N  Average  evaporation  rate  per  unit  area=  M  /  A 


N  Average  2-D  droplet  evaporation  rate  based  on  stream-wise  dimension 

P  Pressure 

Q  Volume  of  drop 

R  Gas  Constant 

Rc  Sessile  drop  radius  of  curvature 

R  r  Coutant  and  Penski  evaporation  parameter  mass^/time 

R  Coutant  and  Penski  evaporation  parameter  zero  convection  velocity 

Rea  Reynolds  Number  =  ud/v 

r  Computer  step  size  parameter  =AX/(AY)2 

Sc  Schmidt  Number  =  v/D 

Shd  Local  Sherwood  Number  =  Nd/cw  D 


Shj 

Shd 

T 
u 
u + 

Ux 


Average  Sherwood  Number  over  drop  diameter  =  N  d/cw  D 
Average  Sherwood  Number  over  circular  drop  area  =  M  d/Acw  D 

Temperature 

Stream-wise  velocity 

Law  of  the  wall  velocity  coordinate  =  u/uT 


Shear  velocity= 


flM 

fUJ 
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X  Computer  variable=urx/(1000vSc) 

x  Coordinate  in  the  flow  direction 

x,  Statistical  analysis  independent  variable 

x0  Boundary  layer  development  length 

Y  Computer  variable=  y+  / 1 0 

y  Coordinate  normal  to  the  surface 

y+  Law  of  the  wall  coordinate  variable  =  uty/v 

y,  Statistical  analysis  dependent  variable 
Greek  Symbols 


p  Contact  angle  of  sessile  drop 

9  Concentration  thickness 

A  Evaporated  vapor  penetration  distance  in  y 

8  Boundary  layer  thickness 

r|  h/d 

v  Air  kinematic  viscosity 

p  Air  density 

<r  Standard  deviation 

<;  y/A 

Subscripts 

A  Ambient  conditions 

b  Slope 

F  Correction  factor 

H  Channel  height 

HD  Chemical  agent  C4H8C1,S 

i,  j  Computer  indices  in  x,  y  direction 

Seg  Segment 

si  Sublayer 

vp  Saturation  vapor  pressure 

vp,Hy  Saturation  vapor  pressure  of  Hydrazine 

y  deviation  from  best  fit  line 
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